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The Equal Filling Approximation, a procedure widely used in mean field calculations to treat 
the dynamics of odd nuclei in a time reversal invariant way, is justified as the consequence of a 
variational principle over an average energy functional. The ideas of Statistical Quantum Mechanics 
are employed in the justification. As an illustration of the method, the ground and lowest lying states 
of some octupole deformed Radium isotopes are computed. 
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I. INTRODUCTION 

The Hartree-Fock-Bogoliubov (HFB) approximation is the cornerstone of the microscopic description of the atomic 
nucleus as it encompasses in the same approximation the concept of mean field orbits needed to understand the 
extra stability associated with magic numbers as well as the concept of pairing correlations needed to understand, 
among others, why the ground state of all even-even nuclei always have the J'^ = O"*" quantum numbers. The HFB 
approximation has widely been used over the last decades with great success both to describe known properties of 
nuclei and to predict the properties of unknown or yet to be experimentally studied nuclei [J, 01 • 

Another nice feature of the mean field approximation is that it allows for the "spontaneous symmetry breaking" 
(SSB) mechanism in which the approximate (mean field) solution of the problem breaks the underlying symmetries 
of the nuclear Hamiltonian. By the SSB mechanism many correlations can be incorporated into the mean field wave 
function while maintaining the simplicity of the mean field description. A minor (mainly practical) drawback of the 
SSB mechanism is that the breaking of symmetries naturally leads to the appearance of "full matrices" in the numerical 
implementation of the method. Those "full matrices" are a direct consequence of the mixing of quantum numbers 
that otherwise could be used to bestow a block structure to the matrices considered. The increase of the effective size 
of the matrices leads to an increase in the number of operation needed to accompHsh the numerical implementation 
of the method and therefore leads to an increase in the computational cost of the problem. Depending on the type of 
problem and the symmetry broken, the computational cost can increase so dramatically as to prevent the large scale 
calculations needed to describe stellar nucleosynthesis or the stability of superheavy nuclei, just to mention a couple 
of physical situations of interest nowadays. A typical case of a computationally demanding application is the study of 
fission barriers allowing for time reversal symmetry breaking, which is characteristic of high spin states or odd mass 
nuclei which have to be treated in the framework of the standard blocking method. In the treatment of odd mass 
systems we have an additional source of computational complexity and it is the self-consistent character of the HFB 
equations, which does not grant that blocking the quasiparticle with lowest excitation energy will yield the lowest 
energy self-consistent solution. As a consequence, several blocking possibilities have to be considered multiplying the 
computing time by the number of the possibilities considered (typically three or four times for each possible Jz and 
parity value). 

On the other hand, the description of odd-A nuclei has started to receive the attention it deserves (three quarters 
of all accessible nuclei are odd-A ones) as the typical odd-even effects, which are not well understood, are intimately 
related to pairing properties and can also serve as more stringent guidelines to the development of new energy 
functional (see for instance 0,0] for recent discussions on this issue). As a consequence, any attempt to reduce the 
computational cost in the theoretical evaluation of odd-A nuclei properties is very important as it will eventually 
allow for a systematic and routinely evaluation of the required properties all over the periodic table. 

A way to reduce the computational cost of mean field calculations when deahng with odd mass nuclei is to try 
to keep time reversal symmetry facilitating in this way the imposition of axial symmetry. To keep time reversal 
symmetry when dealing with odd mass nuclei one is forced to adopt phenomenological approaches in which the 
unpaired nucleon is treated in an equal footing with its time reversed companion. From a practical point of view this 
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phenomenological approach amounts to look at the unpaired nucleon as siting half in a given orbital and the other 
half in the time reversed partner (in the case of preserving spherical symmetry where the orbitals have the 2j + 1 
degeneracy the unpaired nucleon is distributed among all possible angular momentum projections m = —j, . . . ,j with 
equal probability l/(2j + 1)). The above procedure is usually referred in the literature as the Equal (or uniform) 
Filling Approximation (EFA) and has been used quite often in the description of odd nuclei at the mean field level 
and with different interactions - see Refs. [1, 0, 0, S| for recent appHcations of the method. This procedure is used 
because it is considered as an "intuitive" and "reasonable" approach but it is phenomenological in character and lacks 
a solid foundation as there is no product wave function that can reproduce the density matrix and pairing tensor 
of the EFA. The purpose of this paper is to show that the Equal Filling Approximation (EFA) can be described in 
terms of a mixed state (in the sense of quantum statistical mechanics) density operator and the equations to be solved 
are a direct consequence of the variational principle over the energy of such mixed state. As a consequence of this 
microscopic justification it is now possible to introduce numerical procedures like the gradient method to solve the 
EFA equations facilitating enormously the procedure specially in the case of many constraints. Another consequence 
is that now other methods beyond mean field like the calculation of collective masses or the Generator Coordinate 
Method itself can be consistently implemented in the EFA framework. 

Obviously, the EFA is an approximation to the correct treatment of odd-A nuclei in the context of the mean field 
(the blocking procedure, see [2| for a detailed explanation). To discuss the possible differences between them it is 
convenient to look at the odd-A systems as made of an even-even core plus an unpaired nucleon (or quasiparticle) . 
The interaction between the unpaired nucleon and the even-even core will induce polarization effects in the core of 
three types 0, 0, S ED; [lH; namely the mass polarization, the deformation polarization and the spin polarization 
effects. The only one associated to the breaking of time reversal invariance is the spin polarization and it is obvious 
that it is not included in the EFA. The other two effects (mass and shape polarizations) are obviously included in the 
EFA to some extent but it is not easy to estimate the possible impact of spin polarization effects on them. It is clear 
that a deeper understanding of the interrelationship between the three effects is needed and hopefully the present 
justification of the EFA will help to clarify the issue. 



II. THE EQUAL FILLING APPROXIMATION 

In the standard HFB method [2| quasi-particle operators are introduced as Hnear combinations of the creation 
and annihilation single particle operators corresponding to an arbitrarily chosen (usually a harmonic oscillator) basis 

Pfi ~ ^ ^ Umfj.Cjji -f Vfnfj.Cm- (1) 
m 

The HFB ground state wave function is defined by the condition of being the vacuum of all the quasi-particle 
annihilation operators, that is P^\4>) = 0. A more concise definition is given by {(p) = Ilp/^A'I'-') where the index /i 
run over all the quasi-particle annihilation operators that do not annihilate the true vacuum |0). The previous results 
will describe the ground state of an even-even nucleus as it can be shown that a paired HFB wave function is a linear 
combination of product wave functions with an even number of particles. On the other hand, odd-particle systems 
are handled by the "blocked" HFB wave functions 

\4>)^s=P;M (2) 

where //s is any of the quasi-particle indexes compatible with the symmetries of the odd-particle system as, for 
instance, the K quantum number (eigenvalue of Jz) or the parity. As the "blocked" HFB wave function is a product 
of quasi-particle operators. Wick's theorem appHes and the energy is given in the usual way in terms of the "blocked" 
normal and abnormal densities pf'^^) and k^'^^^ 

^(mb) ^ Tr[ip(^«)] + ItiP^^^V^'"'^] - ^Tr[A(^«)K(''^)*] (3) 
The normal and abnormal densities are given by 

pi'k'^ = {ms4ckP+, = {v*v^),,, + (f/fcv,%MB - hvb14; J (4) 

and 
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These two matrices obviously violate time-reversal invariance. As a consequence, the HF field 

''II' - Z^^iq'i'qPqq' W 

qq' 

as well as the pairing field 

qq' 

both violate time-reversal invariance making the numerical calculation much more computationally expensive to carry 
out. A way to preserve time-reversal invariance is to use the Equal Filling Approximation (EFA) that amounts to use 
the "average" density 

and the "average" pairing tensor 

«ffc^^ - (i/*c/^) + \ (c/,,,, v^,^^ - u^,^, + c/,p, v^,-^^ - c/fe.p, v^-^^ ) (9) 

that now preserve time reversal invariance as both expressions involve and average with equal weights of the blocked 
level /is and its time reversed and degenerate partner /Is (see below for higher order degeneracy). Intuitively, the 
above densities should correspond to an occupancy of 1/2 for the states /is and T^s- In the next step of the EFA 
framework, it is assumed without proof that the energy is given by the standard HFB expression but using p^^"^ and 
i^EFA iijg^gafi Qf ^jje corresponding densities, i.e. 

Eefa = Tr[<p^^-4] + ixrp^^^p^^^] - ^TriA^^-^^BFA*] (^q) 

Finally, it is assumed that the U and V amplitudes of the Bogoliubov transformation are given as the solution of the 
standard HFB equation 

/ f^EFA ^EFA \ /UV*\_fUV*\fEO\ 

-A^^^* -h^F^* ) [V U* ) ^ [V U* J [O -E ) ^^^> 

where E are the quasi-particle energies. To our knowledge, the two previous assumptions of the EFA, namely that 
the energy is given by Eq. lfTO|) and that the U and V amplitudes are given by Eq. ([lT|), lacked a foundation and 
were just considered as a plausible quantity (the energy) and equation. Here we will show that both assumptions 
are well founded in terms of standard quantum mechanic procedures and therefore we are giving more credit to the 
approximation. 



A. Justification of the EFA expression for the energy 

In the standard HFB theory the density matrix and pairing tensor are the components of a bipartite generalized 
density matrix 

n4 .^.)4i';;MrMiiVr)-w^w^ m 



(13) 



-K* I- p* J \V U* ) \ Q I ) \ u 
where the generalized quasi-particle density matrix 

1(0I4/3J|0) (01/3^410) ) \0 I 



and the Bogohubov super-matrix 



W^[lll] (14) 
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have been introduced. In the EFA case we can also introduce a generahzed density matrix 

pEFA 
-K^PA* 1-p 



EFA I P^^"^ K^^'^ 
^ = _^EFA* 1 _ ^EFA* (15) 



that again can be written as 

n''^^ = T4^M^^^W^+ (16) 

with 



bEFA _ I Jfj. 



i-U 



(17) 



and the /^^ is given by 



f = } h = Ms or 

■'^ 10 otherwise ^ ^ 



The above result immediately remind us of the Finite Temperature HFB formalism [I4I where the quasi-particle 
density matrix has exactly the same form as above but with the statistical occupancies 



Therefore, the EFA can be viewed as a Finite Temperature HFB formalism with the statistical factors of Eq. (|T8|) . The 
finite temperature formalism is nothing but a quantum mechanics statistical formalism where instead of pure states a 
statistical admixture of them is considered weighted with given probabilities. In the finite temperature formalism the 
probabilities are obtained according to the statistical "ensemble" considered but in the EFA they are just fixed by the 
requirements of the approximation. For this reason we will not use in the following the language of finite temperature 
but instead the one of statistical quantum mechanics. The two relevant concepts in statistical quantum mechanics 
are the one of the "density matrix operator" and the other is the concept of trace. In the present context the trace is 
taken over the whole Fock space in such a way that given a set of quasi-particle creation and annihilation operators 
and and the corresponding vacuum |(/)) (see Eq. Ilj ) we have the following expression for the trace of an 
arbitrary operator 

Tr[6] = ((/-I O 10) + ^ (01 f3^dpl 10) + ^ E (-^1 (^^^-^OplPl 10) . . . (20) 

fj. ' ffJ, 

On the other hand, the density operator V can be chosen in such a way that P|0) = 10) and X*/?!; = Pf.fil'D where 
Pf_i is the probability of the one-quasi-particle excitation 1 0) . In this formalism the statistical mean value of an 
operator is given by 



{0)s = 

with 



= \ ( (01 O 10) + (01 a^Oal 10) + \Y.V,V. (01 /3p/3.6/?t/3t 10) . . . j (21) 



Tr[P] = Z=1 + Y,P^ + Y. P'^P- • • • = 11(1 + p^) (22) 

It is also easy to show that 

{f3pf3a)s = {PjPi)s = (23) 

and also 

(/3+/3,)5 = ^pa-^ - 5p„U; {pppt)s - 5pa (l - = - /,) (24) 

-L ~r Pa \ ^ \ Pa J 
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and therefore we recover the EFA's density matrix of Eq. IjlSp by using the above formaUsm with 

^ \ otherwise 

Here we are impUcitly assuming that the single particle levels are doubly degenerate (Kramer's degeneracy) but in 
those cases where spherical symmetry is preserved in the mean field procedure we will have to populate with the 
same probability all the states with different m = —j, ■ ■ ■ ,j (third component of the angular momentum) for a given 
orbital labeled with the j quantum number. The formalism being developed here apply equally well in this case and 
the reader is referred to Appendix A for technical details in the spherical case. 

Thanks to the existence of a statistical Wick's theorem (see for instance the proof given by Gaudin [l^ or [l3| for 
a more recent account) it is possible to compute any statistical mean value of a product of creation and annihilation 
operators in terms of the corresponding contractions and therefore it is possible to express the statistical mean value 
of the energy {H)s = Tr[/f2?]/Tr[27] by using the standard expression 

{H)s = TrM + iTr[rp] - iTr[A^*] (26) 
where the density matrix and pairing tensor are given by the contractions 

Tr[c+CfcP] Tv[ck-cuV] 

Pkk' = — Kfefc' = -7- — (27) 

Tr[X>] Tv[V] 

Applying this result to the EFA case, we conclude that the energy of Eq. (fTOj) can be written as Eefa = 
Tr[iJ2?^^^]/Ti-[P^^'*]. This result justifies the, otherwise ad-hoc, expression of the EFA energy and gives a physical 
interpretation to it as the statistical mean value of the Hamiltonian taken with the EFA density operator. The EFA 
energy Eefa can also be written in a more transparent way by using Eq. (|2T|) together with Eq. l(25|) as 



Eefa 



\ ((01 H \^) + (01 P^^Hpl^ 10) + (01 ^,HPI^ 10) + (01 fi.sf^sHpljU \^)) (28) 



showing that it is simply an average with equal weights of the energy of the reference even-even wave function |0), the 
energies of one quasi-particle excitations with quantum numbers fiB and Jig and the energy of the two quasi-particle 
excitation with the same quantum numbers. This result, which was to the knowledge of the authors previously 
unknown pgll. is very illustrative of the nature of the EFA as a statistical theory. The same kind of arguments can be 
applied to compute mean values of any kind of operators in the EFA framework. A curious result that can be easily 
derived is that the EFA mean values of any one-body operator, which according to the general result can be written 
as 

(6) EFA = \ ((01 6 10) + (01 a^^Oal^ 10) + (01 a^,64^ |0) + (0| a^^aj^^daj^^al^ |0)) , (29) 
can also be written in a more compact way as 

{d)EFA = I {{<t>\(^^.BOaU \^) + {<f'\<^TZs04s 1'^)) = I {(^\^\^) + (<^l"MB«MB04.aL l<^)) (30) 
This allows us to write the density matrix and pairing tensor as an average over one quasi-particle excitations 

Pkk"^^ = I ((0|a^«c+Cfea+J0) + (0|ap,c+Cfea+J0)) (31) 

and 



""^^ - - (('^l«Mi.Cfe'Cfca+^|0) + (0|a^^Cfe'Cfea+^|0)) (32) 



kk' 
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which is a very intuitive result according to the expressions of Eqs. ([8]) and ([9]). This result, however, 
does by no means imply that the energy, which is the average of a two-body operator, could be written as 

\{{<l>\a^nHa+J4) + {<t>\a-^^Hal^\ 
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B. Variational derivation of the EFA-HFB equation 



Another interesting feature of the results obtained so far is that the variational principle can be applied now to the 
Eefa energy. As it will be shown below the variational principle leads naturally to the HFB equation of the EFA 
mentioned in Eq ifTTj) justifying thereof its use to determine the coefficients of the BogoHubov transformation. This 
result is also advantageous from a practical point of view as the variational origin of the HFB-EFA equations allows 
the use of "gradient-like" methods to solve it and also makes the treatment of constraints much easier. 

The first step in the application of the variational principle is to establish the variational space by defining the most 
general Bogoliubov transformation. This is a common procedure that the interested reader can consult in the standard 
literature [3,[lB| and we will give here only the most relevant formulas just to establish the notation. Given a reference 
HFB wave function \(j)) the most general HFB wave function not orthogonal to it is given by |0(Z)) — exp(iZ)|(/)) 

where Z is an hermitian (to preserve the unitarity of the transformation) one-body operator Z = ^ Z^i^a+aj, 

which is here written in terms of the generahzed quasi-particle operators — {(3i, . . . , (3n , t ■ ■ ^ Pn)^ its hermitian 
conjugate and the bipartite hermitian matrix 

where Z^^ (m, n = 1, . . . , A^) is an hermitian matrix (A^ free parameters, complex Z^^ with m > n plus real Z^^^ ) 
whereas Z^^ is skew-symmetric {N^ — N free parameters, complex with m > n ). The matrix elements of Z^„ 
with m > n plus Z^^^ and those of Z^^ with m > n constitute the complex variational parameters of our model. The 
complex variational parameters will also be denoted by the vector Zp of dimension 2A^ — A (see appendix B for details). 
The coefficients of the BogoHubov transformation of the quasi-particle operators associated to |(/'(Z)) are given by the 
matrix W{Z) = W{0) exp(iZ) which is written in terms of the exponential of Z and the Bogoliubov transformation 
coefiicients 1^(0) of the reference HFB wave function \4>). To determine the dependence of the generalized density 
with the variational parameters we will use Eq. lfT6| to write 



7^^^^^(Z) = WiZ)R^^'^W+{Z) (34) 

where we have kept R^^^ fixed as in Eq. lfT7|). according to its definition. For infinitesimal variational parameters 
(i.e. infinitesimal Z) we have 

n^^^iz) = n^^^iQ) + iw{o)[z, m.^^^]w+{Q) + o{z^) = 7^^■^^(o) + i[z, n^^^] + ©(z^) (35) 

where Z = W(0)ZW^{Q). Now to facilitate the manipulation of different quantities we write the energy as 

Eefa = ^Tia [{H + T)S] (36) 

in terms of the Hamiltonian matrix H = ( j' kinetic energy matrix T — j and the 



matrix '^^^^[^qj^J^I^ ^* p* j derived from the generalized density. The trace is taken over the double 

size space where bipartite matrices are defined. In order to arrive to the above expression we have made use of the 
properties Tr[r2pi] = Tr[rip2] and Tr[A2K*] = Tr[AiK2]* where and stand for the Hartree-Fock and pairing 
fields computed with the density matrix pi and pairing tensor Ki respectively. The two previous relations can be 
written using bipartite matrices as 

Tr2 [{Hi - T)S2] = Tr2 [(H2 - T)S{\ (37) 
This result allows to write the variation of the energy in a more compact way, namely 

SEefa = \ (Tis [[n + T)5S] + Tr2 [[SH - T)S]) = iTr2 [HSS] = ^Tr2 [[7^, + ©(Z^) (38) 

where we have made use of Eq. I|35p and the fact that 6S = STZ. The variational condition SEefa = has to be 
handled with care as not all the parameters of the bipartite matrix Z are variational parameters but as it is shown in 
appendix B the variational condition can be written as [TZ, H] = which is the standard form of the HFB equation. It 
has to be emphasized again that the previous form of the HFB equation in the EFA was just an assumption and now 
we are able to justify it in our framework by simply invoking the variational principle. The HFB equation could be 
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solved using iterative methods as it is customary with the standard HFB equations but we have found more convenient 
to take advantage of the variational origin of the equation in order to use other methods to find its minimum, like the 
gradient method. To implement the gradient method it is more convenient to write the variation of the EFA energy 
in the "quasi-particle" basis as 

SEefa = ^Tr2 [[R, H]Z] + ©(Z^) (39) 

with H =W~^ {0)?iW (0) = i * jjii * ] ■ Using this form and the definitions of appendix B we can finally write 

SEefa — X]p=i ^ i.9E)pZp + O(z^) where {gE)p are the components of the gradient of the energy with respect to 
the variational parameters Zp. In the spirit of the gradient method, by choosing Zp = —rj {g*E)p we make sure we gain 
energy at least to first order in z if the scaling parameter (or step size) 77 is always chosen to be positive. The previous 
election implies that Z — irj^&.W^. The step size t] is estimated in each iteration as to make second order terms 
smaller enough compared with first order ones so that the energy always decreases. Once the Z parameters have been 
determined the wave function W{'L) is computed by evaluating the exponential e*^ by means of a Fade approximation 

to the exponential of the form e^' = Npp{x) / Npp{—x) with Npp{x) = X]fc=o '-fe'^-^'^ (that is, the polynomials in both 

(p) 

the numerator and denominator have the same degree; the coefficients are determined by the standard recurrence 
relation c[f''' = c'^li (2p+i-fc)fc ) ^^^^ has the nice feature of preserving the unitarity of the exponential when used with 
anti-hermitian exponents as it is the case. Usually, a Fade approximation of order p— I suffices and in our numerical 
implementation we have taken e*^ « (1 + 5Z)(1 — ^Z)~^. 



C. Dealing with constraints 

In dealing with constraints, we face the numerical problem of minimizing a given function of the variational param- 
eters EefaC^) subject to some constraints of the kind 

^ = ,. (40) 

where Qi are one-body operators which can be usuafiy written as Qi — J2kk' iQi)kk' ^k^k' (extending the results 
below to the case of operators not commuting with the particle number one is straightforward). Introducing the 

bipartite matrix Qi — (^^^ q* ^ we can write 

* = InQ.S] (41) 
that yields, in analogy to Eq l(38|). to the following expression for the variation of qi 



5q, - ^Tr2 [[7^, Q,]Z] + 0{Z^) = ^Tra [[R, Q,]Z] + 0{I?) (42) 

To consider the constrained minimization procedure we will proceed in the standard way by introducing Lagrange 
multipliers \i and a new functional to be minimized, namely E'^pj^{1) = EefaC^) ~ J2i '^i^ii'^)- The Lagrange 
multipfiers are determined as to make the gradient of i?^^^(Z) orthogonal to the ones of the (?i(Z)'s that will be 
denoted by (ffgj^. Taking into account that the gradient of E'pp^{Z) is given by the vector {gE')p — {9E)p — 
J2j {gqj ) p the orthogonality condition yields 

^^-T.^^^''^^ (43) 
j 

where 

^..-E(5«.);(ff^.)p (44) 
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and 

d^^Y.^9E%{9,,)^ (45) 
p 

Using the explicit form of the corresponding gradients (see appendix B) we can finally write the two quantities above 
as 

% = jTr2[[Q„M][Q„M]] (46) 

and 

- |Tr2[[Q,„M][H,M]] (47) 

It is also convenient to estabHsh a procedure to readjust the constraints as it is customary that in the iterative process 
their values qi sHghtly depart from the desired ones qf\ that is qi = g-^'' +5qi. In the context of the gradient method 
such a procedure is elemental and we only have to replace the chemical potentials \j by \j + 5\j and impose that the 
variation of the values of the constraints 5qi given by Eq. (|42l) yields the desired value. The result is r/JAj = S^^Sqi 
where the matrix elements Sij are the same as in Eq. (|46l) . 



D. Density dependent interactions 

For density dependent interactions, like the Gogny force [15| used in the next section to illustrate the whole 
procedure, we have to define the explicit form of the density dependent part of the interaction for statistical averages. 
It seems natural that, if for a pure state (an HFB mean field wave function in this case) the DD part of the interaction 
is a function of the density of the pure state, then for a statistical average the DD should be the same function but 
of the density of the statistical average; that is, a function of 

p{R) = ^^^1^ = I ((<^i p^R) \<» +T.p^ (<^i MR)Pl \^) + ^Y.p^P'^ (<^i PMmUl \^) ■ ■ • j (48) 

(see Eq. l(22|) for the definition of Z). This prescription has been the one used in previous calculations with the 
Gogny force at finite temperature (ItI. [l8l | as well as by other authors with other density dependent interactions like 
several variants of the Skyrme one [19| . This prescription has the right limit when the probabilities go to zero (pure 
state) and also produces consistent results when the one-quasi-particle energies are computed as partial derivatives 
of the energy with respect to the probabilities: when the above prescription is used the expression for the one-quasi- 
particle energies includes the rearrangement term present in all the HF or HFB calculations with density dependent 
interactions |l5|. It is obvious that for a consistent treatment of the problem, the variation of the energy with respect 
to the variational parameters has to take into account also that the Hamiltonian depends upon them via the DD 
term and the corresponding rearrangement terms have to be considered (see US] for details) . To summarize this 
section, in the EFA case we use the density 

(49) 

for the density dependent part of the Gogny interaction. 



III. RESULTS 



To show an example of the proposed method we have computed the spectrum of several odd-A isotopes of the 
Radium in the range of A between 221 and 231. As it is well known, some isotopes of Ra are known to display 
octupole deformation j2l| in their ground state and therefore, in order to study their spectrum, we will carry out 
calculations constraining the octupole moment to locate the different minima and to check their depths as they 
are relevant for the stability of the configuration against octupole fiuctuations. We will limit the calculation to 
axially symmetric (but refiection symmetry breaking) configurations and therefore each of the blocked levels will be 
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Figure 1: (Color online) Potential energy surfaces as a function of the octupole moment Q3 (in units of b^''^ = lO^fm"^) for the 
odd-mass isotopes of Ra considered and blocking in each of the relevant Jz channels from 1/2 till 11/2. 



characterized (and labeled) by its Jz value (but not parity) . The calculations are performed in the framework of the 
EFA with the finite range and effective interaction of Gogny As it is customary in all the mean field calculations 
with the Gogny force, we have subtracted the kinetic energy of the center of mass motion from the Routhian to be 
minimized in order to ensure that the center of mass is kept at rest. We have also dealt with the exchange Coulomb 
energy in the Slater approximation and neglected the contribution of the Coulomb interaction to the pairing field. 
For the Gogny force we have used the parameter set known as DIS that was adjusted more than twenty years ago 
[13, order to reproduce basic nuclear matter properties and the binding energies of several magic nuclei. The 

HFB wave functions have been expanded in a Harmonic Oscillator (HO) basis containing 14 major shells which is 
enough as to grant convergence in the excitation spectra obtained. 

Due to the self-consistent nature of our procedure it is by no means granted that starting the iterative procedure by 
blocking the quasi-particle of lowest energy the minimization process in going to end up in the lowest energy solution. 
For this reason one has to repeat the minimization process several times using different quasi-particle configurations 
each time (usually ones with the lowest one quasi-particle energy) for the initial blocking. In our case we have 
repeated each calculation three times implying a computational cost 18 times higher (six values of Jz times three 
starting configurations) than the corresponding calculation in an even-even neighbor. By following this procedure we 
can be pretty sure to have reached the lowest energy solution for all values of the octupole moment and mass number. 
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Figure 2: Lowest lying excitation spectrum for the six isotopes of Radium considered. In each panel three spectra are included: 
the one to the left is the experimental one, the other two are theoretical predictions (including octupole deformation effect in 
the middle and not including it, that is Qz = 0, to the right) 
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In Fig. [T] we show the potential energy surfaces (PES) as a function of the octupole moment for the six Radium 
isotopes considered and corresponding to the blocking of the lowest quasi-particles with values ranging from 1/2 
to 11/2. Higher values are not considered here as the corresponding single particle levels lie too far away from 
the Fermi surface as to be relevant for the lowest energy configurations. By looking at the PES for different isotopes 
we learn that the response to octupole deformation strongly depends on the Jz value of the level blocked as it is 
for instance the case in ^^^Ra, where the PES of some levels show a minimum at = whereas others (like the 
ones with Jz = 3/2 and Jz = 7/2) have an octupole deformed minimum at Qs = 3&'^/^. Another interesting fact is 
that breaking refiection symmetry implies energy gains of up to 2 MeV in some cases as compared to the symmetry 
preserving mean field configuration and this amount of energy can not be disregarded in the evaluation of masses. In 
all the cases, the potential energy wells are not very deep indicating the relevance of considering fluctuations on the 
octupole degree of freedom. This fact was already observed in calculations of the same kind and the same interaction 
but for even-even nuclei [13, UHl ^-n^ the conclusion reached there was that a treatment of the octupole fluctuations 
was needed. One of the advantages of the present formulation of the EFA is the fact that the standard methods to 
incorporate correlations (as for instance the collective Schroedinger equation, see Refs [il,!!^] for details) can be now 
generalized to the present case Another improvement to the present treatment is to consider parity projection 
in the manner discussed in [23| for even-even nuclei; work along this lines is in progress and will be reported in the 
near future [23 |. 

Another consequence of the different responses to octupole deformation of the different Jz blocked conflgurations 
is that the spectrum corresponding to the minimum energy is rather different from the one obtained by restricting 
the system to reflection symmetric conflgurations showing the relevance of the octupole degree for freedom for the 
ordering of the spectrum of these and other odd-A nuclei in the region of the Actinide. The different theoretical 
spectra allowing and not allowing octupole deformation are depicted in Fig. [2] along with the experimental data. The 
inclusion of octupole deformation improves the spectrum for several nuclei like ^^^Ra, ^^^Ra and ^^^Ra and makes 
it to look much closer to the experimental one. In fact, the inclusion of octupole deformation on these nuclei allows 
for a correct prediction of the spin of the ground state. For the other cases the inclusion of the octupole degree of 
freedom leaves the spectra more or less unchanged as compared to the Qs = results. On the other hand, it has to 
be kept in mind that when the octupole moment is allowed to take values different from zero, parity mixing is also 
allowed in the wave functions and therefore the different levels lose parity as a quantum number. To restore parity 
symmetry a projection onto good parity is required that would lead to the appearance of two levels with opposite 
parity for each one of the levels breaking the parity symmetry. The energy splitting between the two levels strongly 
depends on the octupole deformation but we can state as a rule of thumb that the energy splitting is going to be rather 
small (a few tens of keV at most) in most of the cases and it will hardly exceed 0.5 MeV. Fortunately, the formalism 
developed in this paper can be extended to the situation of symmetry restoration by means of parity projection and 
the results as well as the whole formalism will be published elsewhere. For the present purposes the only relevant 
information needed is that parity projection will lead to a parity doublet with a not so big energy splitting. A bigger 
splitting could eventually be obtained by fully treating octupole fluctuations as mentioned above and again one of 
the advantages of the present formulation is that the collective masses needed for such a task can be consistently 
evaluated in the present framework. Concerning the comparison with the experiment, we have to keep in mind the 
strong sensitivity of the spectra to tiny details of the underlying single particle states that makes quite diflicult to 
obtain the experimental spectrum in the right order. The accuracy of modern effective interaction only entitle to look 
for an agreement in the number of levels and Jz values in a range of 1 or 1.5 MeV above the ground state but it does 
not entitle whatsoever to sought for an agreement in the ordering of the levels. However, the inclusion of the octupole 
degree of freedom allows for a correct description of the spin of the ground state in flve of the six considered nuclei. 
With this in mind we can conclude that the agreement with experiment is quite good in the whole isotopic chain. 

A procedure that is sometimes used to describe odd-A nuclei is to neglect explicitly all kind of polarization effects 
and treat the quasiparticle excitations in a perturbative fashion 01 • To this end, a reference HFB wave function \fR) 
is computed assuming that is fully paired (that is, is a linear combination of wave functions with even number of 
particles) but the number of particles is constrained to be odd on the average. The wave functions of the ground 
state and excitations of the odd nucleus are built as one-quasiparticle excitations built on top of the reference HFB 
wave function (/3+|(/3fl)). The excitation energies are then given by the corresponding one-quasiparticle energies 

computed as the mean value of the Routhian — {ipu\(3^{H — XN)l3'^\ipR) in an attempt to correct those energies 
perturbatively for the fact that particle number differs from the right value by the quantity N^^j^ = {ipji\l3^iN (3^\(fB) . 
The perturbative correction works well when 7V^^ is small but it is not so reliable when this quantity is large, as the 
chemical potential is usually a few MeV. In Table [J we present the results of such perturbative calculation for the 
nuclei ^^'^Ra and ^^^Ra. The spectrum looks rather similar to the selfconsistent one obtained in the EFA framework 
(see Fig. [2]) but the perturbative one is much more compressed. In this table we also include the values of A'^j', and we 
observe very big values of the order of one in absolute value. These large values together with the neutron chemical 
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2^^Ra 


J. 


(MeV) 




J. 


£m (MeV) 




3/2 


0.000 


0.06 


1/2 


0.000 


-0.09 


5/2 


0.059 


-0.56 


3/2 


0.211 


-0.62 


1/2 


0.216 


0.64 


5/2 


0.410 


-0.81 


7/2 


0.695 


0.94 


7/2 


0.484 


0.89 


1/2 


1.077 


0.93 


1/2 


0.573 


0.89 



Table I: Perturbative results for the nuclei Ra and Ra. The lowest five states in each case have been included. 
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Figure 3: (Color online) Dipole moments as a function of Qa for three representative nuclei and all the Jz values considered 



potential energies (-5.24 MeV for ^^'^Ra and -4.97 MeV for ^^^Ra) make a too big perturbative correction. It has to 
be mentioned that the value of N^\_^ov the ground state is rather small in agreement with the motivation for the 
introduction of \(Pr) given in Ref. m. We can conclude that the polarization effects accounted for by the EFA are 
rather strong and the perturbative treatment, although a reasonable qualitative approximation, is not good at the 
quantitative level. No other mean values and/or physical quantities are considered in this comparison as it would 
imply the evaluation of the perturbative correction due to particle number departures from the physical values (that 
is the evaluation of "chemical potential" like quantities for mean values of arbitrary observables) and this is out of 
the scope of the present work. 

Another relevant physical quantity for octupole deformed nuclei is the intrinsic dipole moment Dq that is directly 
related to the strong E\ transition probabilities observed in these nuclei. It is given as the mean value of the dipole 
operator 

NZ 

Da = e—^{{z)prot - {z)neut) (50) 

in terms of the mean value of the z coordinate for protons and neutrons. The theoretical results for such quantity 
and for each blocked configuration and as a function of Qa are presented in Fig. [3] for three representative nuclei. 
In the first nucleus ^^^Ra the global tendency for Do is to increase with increasing octupole moment for all possible 
configurations with different Jz- On the other hand, for the nucleus ^^^Ra the dipole moment Dq first decreases going 
to negative values and afterward increases to reach positive values (or negative but rather small) as the octupole 
moment increases. Finally, for the nucleus ^^^Ra the dipole moments steadily decrease with increasing octupole 
moment and reaching absolute values greater than in the previous case. The behavior is rather similar to the one 
of the neighboring even-even Radium isotopes as can be observed in Ref. [2^ . This behavior was related in |2^ to 
the increasing occupancy of the neutron jis /2 orbital with increasing number of neutrons and leads to the prediction 
of a minimum in the absolute value of the dipole moment |Do| for the nucleus ^^''Ra. As a consequence, we expect 
substantially lower values of |Do| for the nuclei ^^'^Ra and ^^^Ra as is indeed the case. 

The numerical values of the dipole moments obtained at the minima of the potential energy curves of Fig. [T]are given 
in Table [nl The values are given not only for the theoretical ground state but also for other close lying configurations 
as our theoretical prediction not necessarily coincides with the experimental assignments. Experimental values taken 
from the compilation of [2l| and also from [2S] tell us that for ^^^Ra the J'^ of the ground state is 5/2+ with a value 
of Dq — 0.36 ± 0.10 e fm which is in good agreement with the value of the lowest lying 5/2+ theoretical state. 
For the ^^'^Ra nucleus there are values not only for the ground state but also for some excited ones; the values of Dq 
are 0.129 ± 0.009 e fm, 0.035 ± 0.005 e fm and 0.076 ± 0.004 e fm for the 3/2+(ground state), 5/2+ and 1/2+states, 
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'"'Ra 


"'^''Ra 


""^Ra 


"''Ra 


1 
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0.331 


0.056 


-0.253 


-0.413 


-0.226 


0.000 


3 
2 


0.255 


0.175 


-0.263 


-0.275 


-0.357 


-0.374 


5 
2 


0.427 


0.050 


-0.162 


-0.053 


-0.077 


-0.057 


7 
2 


0.403 


0.254 


-0.098 


-0.471 


-0.415 


-0.279 



Table II: Dipole moments Dq in e fm for the six nuclei considered and obtained for the configuration corresponding to the 
minimum of each blocked configuration with varying Jz value. Only the Jz values up to 7/2 have been considered as this is 
the maximum value of that quantity for all the low lying states. 

respectively. As it can be checked in Table [ll] the agreement between theory and experiment is very satisfactory for 
the three states. In the nucleus ^^^Ra the experimental value is 0.14 ± 0.02 e fm that corresponds in reality to 
the absolute value of that quantity (Dq is extracted from B{E1) values where it enters squared and therefore the sign 
can not be determined in that way). Taking this in account, we can say that there is a reasonable agreements with 
the theoretical prediction which, on the other hand, can be quite strongly affected by fluctuations in the octupole 
degree of freedom Finally, for the ^^''Ra isotope the experimental value is 0.099 ± 0.003 e fm for the 3/2+ ground 
state and this value again agrees reasonably well with the theoretical prediction. The agreement between theory and 
experiment can be considered as quite good specially taking into account the fact that no information on this kind of 
physics was included in the fltting procedure of the force afterwards. 

IV. CONCLUSIONS 

A prescription for the treatment of odd mass nuclei in a time reversal preserving mean field (HFB) framework usually 
known as the Equal Filling Approximation has been justified in terms of standard procedures of quantum statistical 
mechanics. It turns out that the EFA can be described as a mixed state where the blocked one-quasi-particle state 
and its time-reversed counterpart have probability one whereas the others have zero probability. As a consequence, 
the EFA energy is given by an average involving the energy of the underlying even-even system, the energy of the 
blocked one-quasi-particle configurations and the two quasi-particle excitation built out of them. As the energy now 
has a well defined expression in terms of the HFB wave functions it is possible to invoke the variational principle to 
obtain the standard EFA-HFB equation and allowing for the use of more sophisticated numerical techniques, like the 
gradient method, for its numerical solution. The method has been appHed to the study of odd-A Radium isotopes 
as a function of octupole deformation and with the Gogny DIS force and the agreement obtained between theory 
and experiment is quite reasonable. One of the advantages of the present method is the preservation of time-reversal 
symmetry that reduces substantially the computational cost of mean field calculations of odd-mass nuclei. Another 
advantage of the justification obtained in this paper is that the procedure can be extended beyond mean field in a 
consistent way increasing its range of appHcability. 
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Appendix A: HIGHER ORDER DEGENERACY: THE SPHERICAL CASE. 

In the spherical case, the levels to be "blocked" are characterized by angular momentum quantum numbers js^rnB 
with niB = —js, ■ ■ ■ tJb- In the spirit of the EFA each of these levels will be uniformly populated with a fraction 
1 / (2js + 1) of a nucleon what means that in this case the density matrix and pairing tensors in the EFA approximation 
are given by the "average" density 

1 

Pkk'^ = {^*^'^) kk' ^ Yi +T ^ {^k'p.B,]B,niBUl^gj^^rnB ~'^k'nB,3B,mB^ktiBjB,m.B) (Al) 

mB=-jB 
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and the "average" pairing tensor 

1 

'^kk' = {V U )i^f^,+-^ {UkfiB,jB,mB^k'tJ.B,jB,mB ~ ^l"fiBdB,mB^knB,jB,mB) (^^) 

mB = -3B 

The EFA occupancies to be used in the density matrix operator are in this case 

J- ^ I 2iiq^ M = MB,iB,TOB ms = -jB,---,jB ^^^^ 

1 otherwise 

where the index fi has been decomposed in the labels Jb corresponding to the total angular momentum, ms corre- 
sponding to the third component of the angular momentum and finally /is which represents the remaining quantum 
numbers. Once the value of the EFA occupancies are established the formalism developed in the main body of the 
paper can be applied straightforwardly and all the formulas can be used verbatim. 



Appendix B: VARIATIONAL PARAMETERS AND GRADIENTS 

Given a reference HFB wave function \(t>) the most general HFB wave function \(t>{Z)), not orthogonal to it, is given by 
|^(Z)) = exp(iZ)l^) where Z is an hermitian (to preserve the unitarity of the transformation) one-body operator Z = 
5 "^fiv^^^f^v which is written in terms of the generahzed quasi-particle operators = . . . , /3jv, ■, • • • , P^)^ 
its hermitian conjugate and the bipartite hermitian matrix 

that parametrizes the Bogoliubov transformation. Not all the 2 x {2NY parameters of this matrix are independent 
as the matrix Z^^}^ has to be an hermitian matrix (with free parameters, the complex numbers with m > n 
plus real Z^^, i.e. A''^ = 2 x {N{N — l)/2) + N ) whereas Z^° is a complex skew-symmetric matrix (with — N 
free parameters, the complex numbers Z^^ with m > n, i.e. N"^ — N = 2 x {N{N — l)/2) ). As customary we will 
consider Z^^ and Z^^ * as independent parameters instead of the real and imaginary parts of Z^^ and will apply the 
same consideration to Z'^^. As a consequence, the variational parameters of the Bogoliubov transformation are Z^^ 
and Zj^^ with m > n, the real parameters Z^^ and finally Z^^ and Z^^ with m > n. The variational parameters 
can be handled in a compact notation by introducing the vector Zp of dimension 2N'^ — N 

' ^mn rn>n 

ram 

Z^rln m>n (B2) 

^mn rn>n 

Z??* m> n 



As obtained in the body of the paper, the variation of the mean value of an observable can be written as 

5a = ^Tt2[OZ] + 0{Z'^) (B3) 
where O = pR, A] . Taking the most general form of the bipartite matrix 



yQ21 Q22 J 

a little algebra gives da as a function of the variational parameters 



(B4) 



Q [ / ^ \ nm mn) mn ' \ mn nm ) mn ' \ nm mn) mn ~^ \ mn nm) mn / ^ \ mm mm) mm 

m ) 

(B5) 
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Taking now into account the expression of O and the fact that A is the quasi-particle representation of the operator 
given by 



^11 ^20 



(B6) 



where is a skew-symmetric and A^^ is a hermitian matrix if the operator is hermitian (as it should be for any 
observable !) we obtain 



Sa — i ^ Al^„^(fn fm)Zlln — fn — fm)Z^m + C.C. 



(B7) 



This expression can be written in a compact way 6a = J2p^pi9a)p by introducing the vector {ga)p which is the 
gradient of the mean value a with respect to the variational parameters 



*^"m(/"-/m) m>n 

m = m 

-iAl°*{l - fn - fm) m>n 

. i^mn(l - fn - fro) m>n 



(B8) 



Finally, it is important to point out that the requirement {ga)p = is equivalent to [M, A] = 0, a fact that is used in 
the derivation of the EFA-HFB equation. 
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